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ABSTRACT 

We investigate the influence of different analytical parameterizations and fit functions 
for the local star formation rate in AMR simulations of an isolated disk galaxy with 
the Nyx code. Such parameterizations express the star formation efficiency as function 
of the local turbulent Mach number and virial parameter. By employing the method of 
adaptively refined large eddy simulations, we are able to evaluate these physical param¬ 
eters from the numerically unresolved turbulent energy associated with the grid scale. 
We consider both single and multi free-fall variants of star formation laws proposed by 
Padoan k Nordlund, Hennebelle k Chabrier, and Krumholz k McKee, summarised 
and tested recently with numerical simulations by Federrath k Klessen. We find that 
the global star formation rate and the relation between the local star formation rate 
and the gas column density is reproduced in agreement with observational constraints 
by all multi free-fall models of star formation. Some models with obsolete calibration 
or a single free-fall time scale, however, result in an overly clumpy disk that does not 
resemble the structure of observed spirals. 
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1 INTRODUCTION 


Numerical simulations of disk galaxies cannot fully resolve 
processes in the interstellar medium (ISM) such as star for¬ 
mation and feedback from supernova (SN) explosions. For 
this reason, effective models for these processes are used in 
large-scale simulations, which are called sub-resolution or 
sub-grid scale (SGS) models. The most commonly applied 
model of star formation is based on the empirical Kennicutt- 


Schmidt (KS) relation 

(e.g. |Schmidt||1959| Kennicutt 

1998 

Gao k Solomon 2004; 

Daddi et al. 2010; Genzel et al. 

2010 

Lada et al.|2010| |Bigie 

et al. 2011|; Kennicutt k Evans 

2012 

Schruba|2013 Roychowdhury et al.|2015). Formulated as an 


SGS model, this relation translates to 


Tff 


( 1 ) 


where p can be either the total gas density (e.g. Renaud et al. 


2013 Agertz et al. 2013| Lagos et al.|2015 Tasker et al. 2015 [ 
or some fraction of that density (for example, the molecular 
hydrog en d ensity as in Gnedin et al. (2009), Dobbs k Pringle 


(2013), or Agertz k Kravtsov (2015)), Tff ~ (Gp) i/z is 


the free-fall time scale associated with p, and eff the so- 
called star formation efficiency. The coefficient eff is usually 


* E-mail: hbraun@astro.physik.uni-goettingen.de 
f E-mail: wolfram.schmidt@uni-hamburg.de 


assumed to be a constant, i. e. a freely tunable numerical 
parameter, in disk galaxy simulations. This constant is cho¬ 
sen such that the global star formation rate matches ob¬ 
servational values. However, numerical studies and observa¬ 
tions of star-forming regions in galaxies suggest significant 


local variations of the star formation efficiency (Evans et al. 
2009 1 [Onodera et al.||20T0| |Murray||2011| |Federrath||2013| 
Sali m et al.|2015 ), which may arise from different evolution¬ 
ary stages of local star-forming environments (Kruijssen k 
Longmore||2014 ). 


Also from the theoretical point of view, the star forma¬ 
tion efficiency should change with the local conditions in the 
ISM because stars form from gravitationally unstable den¬ 
sity enhancements seeded by turbulence in molecular clouds 
(e.g. Mac Low k Klessen 2004| McKee k Ostriker 2007 


Hennebelle k Falgarone 2012). Although the dynamics of 


star-forming clouds appears to be extremely complex (com¬ 
pressible turbulence, self-gravity, magnetic fields, radiative 
cooling and heating), attempts have been made to capture 
the essence of the star formation process in relatively sim¬ 
ple analytical relationships between the star formation rate 
or efficiency and some fundamental parameters. In particu¬ 
lar, star formation efficiencies were parameterized in terms 
of the turbulent Mach number (corresponding to the ratio 
of turbulent and thermal energies) and the viral parameter 
(the ratio of turbulent and gravitational energies). An rep- 
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2 H. Braun and W. Schmidt 


resentative example is the relation proposed by |Padoan~fc| 


Nordlund (2011) (a complete overview will be given in Sec¬ 
tion^, which is based on the idea that the typical length 
scale of dense clumps is given by the thickness of shock lay¬ 
ers for a particular Mach number of turbulence and that 
the clumps become unstable if they exceed their Jeans mass 
(this is where the virial parameter comes in). The distribu¬ 
tion of density enhancements is determined by a log-normal 
density distribution, which again depends on the turbulent 
Mach number ( |Federrath &; Banerjee 2015). Braun et al. 


( 2014| [hereafter B S14| exploited this relation in order to dy¬ 


namically compute the efficiency parameter in simulations 
of an isolated disk galaxy, rather than assuming some spe¬ 
cific value. It turned out that this kind of SGS model is 
able to predict typical star formation efficiencies in galaxies 
(e.g. Gao &; Solomon||2004| Daddi et al.||2010 Genzel et al. 


|2010 Murray 2011) without making use of the observed rela¬ 
tions for the star formation rate as input. However, an open 
question is whether the properties of simulated disk galax¬ 
ies are sensitive to the various assumptions that go into the 
star formation model. If that were the case, it should be 
easily possible to discriminate existing models, which make 
entirely different basic assumptions. On the other hand, it 
is also possible that different models produced comparable 
effects under conditions that favour star formation. 

Another essential component of a SGS model for real¬ 
istic galaxy simulations is stellar feedback. There is a great 
variety of different approaches and extensive studies of the 
influence of feedback models on the properties of simulated 


galaxies (e.g. Stinson et al. 2006 

2013 Wise et al. 2012; 

Agertz et al. 2013 Hopkins et al. 

2014). Here we focus on 


the influence of the star formation model on global prop¬ 
erties of simulated isolated disk galaxies for one particular 
feedback model, namely the mixed thermal and turbulent 
feedback introduced in lBS141 It was demonstrated that this 
model produces reasonable results and therefore can be con¬ 
sidered as a reliable reference point for further studies. 

This article is structured as follows. After outlining the 
applied numerical methods and summarising the star forma¬ 
tion models implemented in our code (Section|2j, we present 
the results from our simulations in Section [3] In particular, 
we consider the global star formation rate, the structure of 
the gaseous disk, the relations between the local star forma¬ 
tion rate and surface densities, and the role of turbulence 
for the different models. The results of our study are sum¬ 
marised and discussed in Section [4] 


2 NUMERICAL METHODS AND MODELS 


The simulations presented in this study were carried out 


using the cosmological hydrodynamics code Nyx (Almgren 


et al. 2013). Nyx solves the Euler equations on an adap¬ 


tively refined grid using the piecewise parabolic method. 
Collisionless components like stellar particles are taken into 
account as N-body system via a particle-mesh approach. In 
our simulations, we included the SGS model introduced by 


Schmidt & Federrath (2011) and Schmidt et al. (2014) to 


model turbulent motions below the resolution scale. Non- 
adiabatic sub-resolution processes are modelled using the 
MIST (Multi-phase Interstellar medium with Star forma¬ 
tion and Turbulence) model, which was proposed in Braun 


& Schmidt (2012 hereafter BS12) and BS14 The employed 


physics encompasses: 


• Self gravity from gas and stars. 

• Gas dynamics using piecewise parabolic method for the 
resolved motions and the SGS model for unresolved turbu¬ 
lence. 

• Radiative cooling from gas, metals, and dust. 

• Multi-phase ISM consisting of a diffuse warm, a clumpy 
cold, and a hot phase for SN ejecta. 

• Star formation rate depending on the locally inferred 
molecular fraction and the thermal and turbulent state of 
the gas. 

• Stellar feedback in the form of a combination of thermal 
and turbulent feedback from SN and thermal feedback from 
Lyman continuum heating, both depending on the age of 
the stars. 

• Metal enrichment due to SN. 


As demonstrated in [BS14[ simulations featuring MIST in 
combination with a SGS turbulence model are able to si¬ 
multaneously reproduce several observed properties of star 
formation in a plausible way. Owing to the combined effect 
of turbulent plus thermal SN-feedback and a non-cooling, 
but decaying hot gas-phase of SN-ejecta the so-called over¬ 
cooling problem is avoided. The star formation relations of 
KS-type following from the simulation data are a genuine 
result of the modelling, as they were not a priori imposed in 
the model. Moreover, MIST allows us to probe the behaviour 
of different theoretical approaches to the calculation of the 
star formation rate in the cold molecular gas by incorporat¬ 
ing the corresponding parameterizations in our simulations. 


2.1 Multi-phase and star formation models 


MIST separates the gas content of a numerical resolution 
element into two fractions - or phases: a clumpy cold phase 
and a diffuse warm phase, within which the cold clumps are 
embedded. The fractional densities of the cold and warm 
phases are determined by partial differential equations with 
source terms related to various exchange processes between 
the phases and star formation, which converts cold gas into 


stellar mass (see Figures 1 and 2 of BS12 for an overview). 
While the warm phase is subject to radiative cooling and 
feeds the clumpy cold phase through this channel, the cold 
phase itself is kept at a constant temperature of T c = 50 K. 
The heating of cold gas is modelled as mass transfer into 
the warm phase, which has a variable temperature following 
form the hydrodynamical energy equation. Under adiabatic 
conditions, the total energy of the two phases is conserved. 
Moreover, shielded molecular gas may reside inside the cold 
clumps. The mass fraction /h 2 of molecular gas is estimated 
using a Stromgren-like approach. 

Only this molecular fraction /h 2 Pc of the cold gas with 
fractional density p c is allowed to form stars at a rate p s 
given by 

_ €core€mod/H 2 Pc 


T c ,ff 


7c,ff is the free-fall time scale within the cold clumps 


7c, ff = 


37T 


32Gp c ,p 


( 2 ) 


(3) 
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How the star formation law affects galactic disk structure 3 


with the gravitational constant G and the average density 
inside the cold phase p c , P a, which differs from the fractional 
density p c . To calculate the average densities, an equilibrium 
between the effective (thermal plus turbulent) pressures of 
the cold and warm phases is assumed at the typical clump 
scale £ c (for a thorough description see BS12[ BS14). The co¬ 
efficient emod in equation |2]) specifies the efficiency of form¬ 
ing gravitationally bound cores with respect to the free-fall 
time scale r c ,ff. The core efficiency parameter e cor e = 0.4 
takes into account that only a fraction of the mass in the 
bound cores eventually ends up in stars. The remaining frac¬ 
tion (1 — e C ore) of the prestellar material is assumed to be 
re-injected into the ISM due to protostellar jets, winds, and 
outflows. 

No explicit density threshold is imposed on the total 
gas density p or the average density of the cold phase, p c ,pa- 
The only requirement for star formation is the existence of 
shielded molecular gas inside the cold clumps, which de¬ 
pends on p c ,pa, £ c , and the metallicity Z. 

A unique feature of MIST is that e mo d is not assumed 
to be a constant coefficient, but is treated as a variable of 
the model, which is computed from the current local state 
of the gas in the cold clumps. This is achieved by mod¬ 
elling the substructure of the ISM and turbulence below the 
resolution scale on the basis of the numerically unresolved 
turbulent energy i^sGS, which determines key parameters in 
the expressions for e mo d following from analytical theories of 
star formation. 

In |BS12| an d |BS14[ the model of |Padoan &; Nordlund 


(2011 hereafter PN11| in the parameterization of |Feder- 


rath &; Klessen (2012 in the following FK12) were used. 


In this study, we explore the effects of different choices for 
e m od based on the analytical models proposed by |PN11[ 


Krumholz & McKee (2005 hereafter KM05), and Hen- 


nebelle & Chabrier (2008 2009 2011 2013 hereafter HC08- 


13). While these models result from theoretical reasoning, a 
simple star formation law was put forward by |Padoan et al.| 
(2012 hereafter PHN12) by fitting data from simulations of 


driven isothermal turbulence in self-gravitating gas. 

In the context of the analytical models, e mo d can be 
understood as an x-weighted integral of the probability 
density function pdf (a?) over the logarithmic overdensity 
ratio x — log(p/p c , P a) above a certain critical threshold 
x cr = log(p cr /pc,pa), which separates gravitationally bound 
from unbound structures in the turbulent gas. To account 
for a collapse time scale that varies with the density fluctua¬ 
tion x , multi free-fall (mff) models incorporate the weighting 
factor /ff = T c ,ff/Tff(p) = exp * 1 / 2 (x) into the integral. In gen¬ 
eral, e mo d can be written as 


COG 

Cmod = /cal / fs(x) exp(x) pdf(x) dx. 

J x cr 


( 4 ) 


The arbitrary prefactor / ca l (1 /<pt in FK12) has to be cali¬ 
brated by fitting the star formation efficiency to numerical 
data from simulations of self-gravitating turbulence. 

Although pdf(x) tends to develop a power-law tail 
in strongly self-gravitating gas (e.g. Girichidis et al. 2014 


Kainulai nen et al.|2009 ), the analytical models generally as¬ 


sume a log-normal shaped PDF (the reasons of which are 
discussed in Federrath & Klessen| ( |2013| ), which is charac¬ 


teristic for supersonic isothermal turbulence: 

1 ( (x — Xo ) 2 


pdf©) — 




exp 


2 cr.2 


where 


( 5 ) 


at = log (l + b 2 M 2 ) , (6) 

determines the width of the distribution and xo = —(t 2 /2 is 
the logarithmic mean. Apart from the compressive factoiQ 


J 1 if local SNe feedback is active 
1 0.5 in quiescent regions 


( 7 ) 


which was introduced by Federrath et al. (2010), a x depends 
on the root-mean-square Mach number A4 of turbulent mo¬ 
tions in the cold clumps. We estimate M in our simulations 
by extrapolating the turbulent SGS energy from the grid 
scale to the clump scale £ c by assuming a power law for the 
scale dependence of turbulent velocity fluctuations: 


/ 2W c/A)** 
V 7(7 - l)e c 


with the specific thermal energy of the cold phase e c = 
e(T c = 50 K), the adiabatic index of the equation of state of 
the gas 7 = 5/3, the specific kinetic energy Ksgs of motions 
below resolution scale A, and the turbulent velocity scaling 
exponent rj w — 1 / 3 . 

Basically, the differences between the analytical models 
amount to different definitions of x cr and /ff. The critical 
overdensity x cr is in all cases a function of M and the virial 
parameter <a v ir of the cold clumps, where 


looses (4/A) 2,) 

7T G£%Pc,p a 


(9) 


In addition to the models in their originally published 
form, we also consider the single and multi free-fall formu¬ 
lations with the best-fit calibrations as provided by |FK12| in 
their tables 1 and 3. A comprehensive fist of all the models 
used in this study is given in Table [l] The main ideas behind 
these models are briefly summarised in the following. For a 
more detailed account see |FK12| and the original papersj^] 


2.1.1 PN family 

|PN11 [ assume that the size of a critical Bonnor-Ebert-sphere 

is equal to the typical thickness of a shocked layer, which 
is computed from the combination of the isothermal shock 
jump conditions and a turbulent velocity scaling relation 
u'(£) oc £ ??c with 7] c = 0.5. This implies 

exp(xcr) = 0.0067 9~ 2 a v i r A4 2 . (10) 


1 This definition assumes that turbulence injection driven due to 
SNe is dominated by compressive/dilatational forces. Otherwise 
a mixture of compressive and solenoidal driving forces resulting 
from shearing motions in the galactic disk, local gravitational 
collapse, or motions triggered by thermal instabilities of the multi¬ 
phase ISM is assumed. 

2 In the following we refer to different subsets out of our simula¬ 
tion suite as families. PN, PN-FK, and PN-FK-mff belong to the 
PN-family, HC, HC-FK, and HC-FK-mff to the HC-family, and 
KM, KM-FK, and KM-FK-mff accordingly to the KM-family. 
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4 H. Braun and W. Schmidt 


Table 1. Overview of star formation models implemented into MIST. For a brief description of the models and definitions of the various 
coefficients and parameters, see Sections [XOpT3l 


ID critical overdensity ratio ff-time factor SF efficiency coefficients references 

exp(> cr ) /ff e mod 


models and coefficients as in original publication (see listed reference - if they did not suggest a value, we chose f ca l = 1) 


PN 

(0.0067 )6- 2 a vir M 2 



exp 1 / 2 (x cr ) 

/cal X a 2 

e = 0.5, /cal = 1.0 

PN11 



HC 

( 71 ’ 2 /5)2/ cu t°Nir (y c ut M 

- 2 + 1 / 3 ) 

exp 1 / 2 (x) 

/cal X a 3 

ycut — 0.1, fe a 1 = 1.0 

HC08-13 



KM 

(7T 2 /5)/ 2 a v irA4 2 



1 

/cal X CL\ 

<t>x = 1.12, /cal = 1.92 

KM05 



PHN 

n.a. 



n.a. 

exp [<p7r 2 a v i r ] 

<p = —0.003 

PHN 12 



calibrated by FK12 (see their table 1 for analytic forms and their table 3 (HD-fit) for the coefficients) 

PN-FK 

(0.0067 )0- 2 a vir M 2 



exp 1 / 2 (a? cr ) 

/cal X (22 

e = 1.0, /cal = 1.5 

PN11 FK12 

HC-FK 

( 7r2 /5)2/ cu t°Nir (y cu t M 

- 2 + 1 / 3 ) 

exp 1 / 2 (x) 

/cal x a 3 

ycut — 1-3, fea 1 — 0.24 

HC08-13 

i FK12 

KM-FK 

(7T 2 /5)/ 2 a v irAf 2 



1 

/cal X Off 

4>x = 0.12, /cal = 3.0 

KM 05 "I 

Tvl2 


multi free-fall (mff), calibrated by 

FK12 

(see 

their table 1 for analytic forms and table 3 (HD-fit) for the coefficients) 


PN-FK-mff 

(O.OO67)0- 2 a vir .M 2 



exp 1 / 2 (x) 

/cal x a 3 

0 = 1.0, / ca i = 0.49 

PN11 FK12 


HC-FK-mff 

(7r 2 /5)y~ u 2 t a vir M~ 2 



exp 1 / 2 (a?) 

/cal x a 3 

ycut — 1.1 5 /cal =0.21 

HC 08-13 

i FK 

12 

KM-FK-mff 

(7T 2 /5)0 2 a v i r M 2 



exp 1 / 2 (x) 

/cal X <23 

<t>x = 0.19, /cal = 0.49 

KM05 I 

r K12 



Definitions: 


“l = 0.5 (1 + erf[(cr 2 - 2x cr )/(8<7 2 ) 1/2 ]); 

<22 = 0.5 (l + erf [(a 2 — 2x cr )/(8<7 2 ) 1 / 2 ]) exp (x cr /2); 
<23 = 0.5 (l + erf[(<r 2 x C v) / (2<r 2 ) 1/2 j) exp( 3 <r 2 / 8 ) 


In the MIST framework, 6£ c is the injection scale of super¬ 
sonic turbulence in a cloud. In the case of our fiducial model 
PN-FK we simply set 0 = 1.0, while 6 — 0.35 in |PNll| and 
0 = 0.65 in table 3 of lFK121 We chose the intermediate value 
6 = 0.5 for our PN run to explore the influence of this co¬ 
efficient. Both the original [P N11 1 model and PN-FK assume 
a constant free-fall time factor /ff = exp 1 / 2 (x cr ), i.e. the ef¬ 
fective free-fall time scale is given by the critical overdensity 
Per rather than the mean density p c ,pa of the cold phase. 


tuting the definitions of M and <a v ir yields 


2 2 

/ \ 7T G!vir , 7T C^vir 

exp(x cr ) = —^-7-7o + 777-• 

W2 15ycut 


%c utM 2 


( 12 ) 


This expression is used for the models HC and HC-FK. Fol¬ 
lowing |FK12[ the second term is dropped for the HC-FK-mff 
model. For the HC family, ’mff ’ is actually a misnomer be¬ 
cause all models in this family feature the multi free-fall time 
factor /ff = exp -1 / 2 ©). 


2.1.2 HC family 


The model suggested by |HC08-13| is based on the assumption 
that bound substructures of a turbulent gas cloud collapse 
on time scales determined by the mean densities of the sub¬ 
structures rather than the whole cloud. The star formation 
efficiency is then obtained by integrating the mass spectrum 
of substructures weighted by their free-fall time factors. This 
integral is truncated at a certain mass fraction y cu t to avoid 
overly large structures. While |HC08M3| suggested y cu t ~ 0.1, 
|FK12| concluded from their fits to simulation data that y CVL t 
slightly above unity is favoured. 

Assuming that additional non-thermal support against grav¬ 
ity must increase the stability range (a notion that was ques¬ 
tioned by Schmidt et al. 2013), the critical overdensity is de¬ 
fined by requiring that the turbulent Jeans length Aj,t(p cr ) 


Aj,t(pcr) — 


71-7(7 - l) e c + 27 rAj,t(p C r)ifsGS (^) 
G p C r 


277^ 


(n) 

equals the size y c ut£c of the largest collapsing substructures. 
This condition leads to a quadratic equation in £ c . Substi- 


2.1.3 KM family 

|KM05[ argue that the critical overdensity follows from the 
condition that the sonic length is close to the Jeans length. 
From this condition it follows that 

exp(x cr ) = (7r 2 /5)0^O!virA4 2 . (13) 

Apart from the prefactor, the resulting expression for 
exp(xcr) is the same as for the models in the PN family. 
The fudge factor <f x should be of the order unity. |KM05| ob- 
tain the best fit for <f x — 1.12. However, |FK12| determined 
smaller values (cp x ~ 0.1... 0.2). 

The main difference between the original PN and KM mod¬ 
els is that there is no free-fall time factor at all in the KM 
models (/ff = 1 in the case of KM and KM-FK). This means 
that |KM05| assume that the free-fall time scale is just given 
by the mean density of the cold gas (see equation [3] for r c ,ff), 
while PN and PN-FK set the effective free-fall time scale in 
the star formation law equal to Tff(y5 cr ) (in this case, r c ,ff 
cancels out from equation © for the star formation rate). 
The mff variants of the KM and PN models, on the other 
hand, are formally the same. 
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2.2 Simulations 


For each of the ten star formation models in Table [l] a simu¬ 
lation run of an idealised isolated disk galaxy was performed 
with Nyx (Almgren et al. 2013) using the methodology from 
(BS14). As initial conditions we used data obtained from 


the ’ref’ of BS14 run with their fiducial model (in this pa¬ 


per referred to as PN-FK) at a simulated age of roughly 
1 Gyr. Starting the simulations for our comparison from a 
pre-evolved galaxy has several benefits. Employing an aged 
galaxy as initial conditions minimises transient effects of the 
rapid initial growth of the stellar mass and the metallicity 
within the galactic disk. Since we begin with a statistically 
stationary star forming configuration of the galactic disk, in 
which star formation and stellar feedback effects are nearly 
balanced, it is easy to detect the changes of the disk config¬ 
uration after switching to a different star formation model. 

The ’ref’ run of IBS14I simulated the evolution of an 
idealised isolated disk galaxy residing in a box of (0.5 Mpc) 3 
size with the PN-FK model. With a root grid of 256 3 and 
6 levels of adaptive mesh refinement an effective resolution 
of 30 pc was achieved. The galaxy in ’ref’ was initialised as 
an adiabatically stable, rotating gas disk of 10 10 M© mass 


using the potential-method described by Wang et al. (2010). 


The static potential of an NFW-shaped dark matter halo of 
10 12 M© is added to the self-gravity of the gaseous disk. 
There was no initial stellar component. Radiative cooling 
causes the initially thick and hot disk to rapidly collapse 
into a fragmenting thin disk. 

After 1 Gyr of evolution, about 30 per cent of the initial 
total gas mass has been converted into stars. While the ma¬ 
jority of the long-lived stars are rather smoothly distributed 
over the inner disk, many young stars reside in stellar clus¬ 
ters with a typical lifetime of a few 100 Myr. The gas content 
of the inner, star-forming part of the disk is below 50 per 
cent. At this stage, the simulated galaxy roughly resembles 
a disk galaxy in a quiescent star forming state with a global 
star formation rate of ~ 1 M©yr -1 governed by the self¬ 
regulation between star formation and stellar feedback(^]We 
use this configuration as initial condition for our comparison 
runs by following the disk evolution over additional 0.4 Gyr 
with different star formation models. The time period of 
0.4 Gyr corresponds to one orbital revolution at 10 kpc dis¬ 
tance from the galactic centre or roughly ten times the max¬ 
imum age of a stellar particle producing feedback due to SNe 
of type II. 


3 RESULTS 

3.1 Global star formation rate 

The star formation rate integrated over the whole disk is 
plotted as function of time in Fig. [I] for all runs. After a 
short transitional phase, in which the disk configuration ad¬ 
justs itself to the new conditions set by the star formation 
model, the global star formation rate Msf settles on aver¬ 
age between 2 M©yr -1 and 3 M©yr -1 in the majority of the 

3 Our configuration is not comparable to a star-burst galaxy at 
high redshift such as the ’HighZ’ runs in [H opkins et aT| (2011). 
After 1 Gyr, the conditions in the inner disk of our simulation 
roughly correspond to the ’Sbc’ runs in that paper. 



Figure 1. Top panel: Global star formation rate Msf versus 
simulation time for runs using the different star formation models 
listed in Table ^ The evolution of Msf in runs with PN11 family 


models is drawn in black, with |HC08-131 family models in green, 
with KM05 -family models in red, and with the |PHN12| model 
in blue. A model in its original version is indicated by a dot- 
dashed line, its |FK12f calibrated version by a solid line, and its 
|FK12f calibrated multi free-fall version by a dashed line. Bottom 
panel: Same as above, but Msf is smoothed on a time scale of 
80 Myr using a moving average to filter out short term variations. 


runs (PN-family, HC-family, and KM-FK-mff). In the bot¬ 
tom panel of Fig. [lj short term variations are smoothed out 
by applying a moving average with width of ~80 Myr. These 
variations, which can be as large as 1 M©yr _1 in amplitude, 
are a consequence of the life-cycle of individual star-forming 
regions in the disks. 

As explained in |BS14[ the associated time scale of 
around 10 to 30 Myr arises from delayed SNe feedback. This 
time scale is consistent with observational estimates of the 


life-time of star-forming gas clouds (e.g. Blitz et al. 2007 


McKee &; Qstriker||2007 Miura et al.||2012 ) and is also re¬ 


produced by numerical simulations (e.g. Hopkins et al.|2012 


Tasker et al.||2015| Dobbs &; Pringle||2013 ). The smoothed 

time evolution of Msf is similar for most runs, including the 
PN-FK run, which merely continues the ’ref’ run of |BS14| 
However, the runs KM, KM-FK, and PHN exhibit de¬ 
viating behaviour. These three runs have in common that 
Msf drops well below 1 M©yr _1 shortly after the start of 
the simulation and then gradually grows until the simulation 
is stopped. Only the PHN run with the simple star forma¬ 
tion model seems to enter a self-regulated regime toward the 
end of the simulation, where Msf approaches a similar level 
as in the other simulations. While the KM-FK model leads 
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Figure 2. Comparison of area-weighted 2D histograms of the star 
formation efficiency e mo d in the cold phase vs. the SGS turbulence 
energy Ksgs at an age of ~ 1.4 Gyr. The area is logarithmically 
colour-coded. 


to an overproduction of stars compared to the continuation 
of the ’ref’ run with the PN-FK model, the original KM 
model falls much below the expected star formation rate. 
The deviations observed for the KM approach are remedied 
with the multi free-fall modification of the star formation 
law. Otherwise, the global star formation rate appears to 
be quite robust with respect to the parameterization of star 
formation. 


3.2 Turbulence 


Figure [2] reveals that star formation occurs for all included 
models preferentially around the sweet spot at specific SGS 
energies of Ksgs ~ 100 km 2 s _1 . This corresponds to a pre¬ 
ferred velocity dispersion of cr u ~ 10 km s -1 on the grid res¬ 
olution scale A « 30 pc, implying turbulent Mach numbers 
around 10 in the cold clouds on the length scale CQ These 
values are comparable to observational findings (e.g. | Leroy 
et al. 2008 Shetty et al. 2012 Stilp et al. 2013). The intensity 


of small-scale turbulence in actively star-forming regions is 
mainly determined by stellar feedback, the exchange of mass 
between warm and cold phases, and the assumed effective 
pressure balance. Figure [2] shows that the models in the PN- 
and HC-families as well as the KM-FK-mff model predict 
roughly comparable star formation efficiencies e mo d ~ 0.1 


4 Note that A4 C is different from the Mach number cr u /c c scaled 
down from A to £ c because the latter encompasses both the warm 
and cold phase. 


(corresponding to star formation time scales of < 100 Myr) 
for the relevant SGS energies around 100 km 2 s -2 , although 
their dependence on Ksgs is quite different. The highest ef¬ 
ficiencies are obtained for the HC-FK-mff model, while PN 
tends to produce relatively low efficiencies. This is also re¬ 
flected by the Esf-Sh 2 -relations discussed in section 3.4 


Moreover, one can clearly see that the efficiencies follow¬ 
ing from the KM, KM-FK, and PHN models are far below 
the plausible range of values in turbulent star-forming re¬ 
gions. For these models, star formation exclusively occurs in 
regions of exceptionally high density, where efficiencies that 
are systematically too low are compensated by much shorter 
free-fall time scales. 

The star formation efficiencies obtained with the differ¬ 
ent models can be understood as follows: 


3.2.1 PN family 


All PN-family models produce a similar e mo d-7FsGS rela¬ 
tion. However, the efficiencies are systematically lower for 
the original PN model and the implied star formation time 
scale is about 0.75 Gyr, much longer than ~ 100 Myr for 
PN-FK and PN-FK-mff. Due to the lower level of e mo d, the 
disk tends to be clumpier. This is mainly a consequence of 
the chosen model coefficient 6 — 0.5 compared to 0 — 1 for 
the FK-calibrated variants of the PN-family (see Table [TJ. 
This behaviour supports our assumption that i c should be 
close to the integral scale of turbulence in cold gas clouds, 
but the strong sensitivity on the parameter 0 makes PN-type 
models difficult to control. 

The declining trend of e mo d with increasing SGS turbu¬ 
lence energy is due to exp(x cr ) oc M 2 oc Ksgs (see equa¬ 
tion 10), i.e. stronger turbulence tends to produce thinner 


shock layers, requiring higher densities to reach the criti¬ 
cal mass. This effect is alleviated, however, by the implicit 
dependence of p c , P a (the mean density of the cold phase) 
on M. This dependence is a consequence of the effective 
pressure equilibrium between the phases. Since mostly loga¬ 
rithmic density fluctuations x close to x cr contribute to the 
integral over the exponential tail of pdf(x) for x > x cr , the 
weighting factor (i.e. single vs. multi free-fall) has only a 
minor impact. 


3.2.2 HC family 

All HC-family models are actually multi free-fall models and 
yield realistic e mo d ~ 0.1. As a consequence, the resulting 
star formation rate is comparable to the PN-FK case, al¬ 
though the trend of e mo d with Ksgs is opposit e to the PN 
family because exp(x cr ) oc 2 (cf. equation 


12 


and 


10 ). 


Omitting the second term in the expression for x cr in the 
case of HC-FK-mff (see Table [TJ significantly affects e mo d 
only for values of Ksgs that are generally lower than in 
star-forming regions. 

Although the values chosen for the cutoff mass y cu t dif¬ 
fer by a factor 10 in the HC and HC-FK runs, we do not see 
a significant impact on e mo d- However, this does not neces¬ 
sarily imply that y cu t is unimportant. Probably the influence 
of 2/cut is compensated by the calibration factor / ca 1 « 0.2 
in the case of HC-FK compared to / ca 1 = 1 for HC. It is 
therefore likely that a careful calibration is also required for 
HC-type models. 
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3.2.3 KM family 

KM-FK-mff differs from PN-FK-mff only by a constant co¬ 
efficient, which explains the striking similarity of the top 
right and the bottom right plot in Figure [2] The functional 
dependence on the Mach number, exp(x cr ) oc Ai 2 , is the 
same for all KM and PN models. However, both models, 
KM and KM-FK lack a free-fall weighting factor in the in¬ 
tegrand: fff = 1, independent of x cr (see second column in 
Table [I]). For this reason, e m od rapidly drops with increas¬ 
ing Ksgs oc Ai 2 , which strongly suppresses star formation 
as soon as stellar feedback kicks in (about 4 Myr after cre¬ 
ation of the first stars). The suppressed star formation in 
turn prevents feedback from becoming sufficiently strong to 
expand and destroy a region with active/recent star forma¬ 
tion. This results gradually growing global star formation 
rates in the KM and KM-FK simulations, without reaching 
statistically stable disk configurations. 

3.2.4 PHN-model 

In contrast to the analytical models, |PHN12| assume a strong 
increase of e mo d with the cold-phase density through the ex¬ 
ponential dependence on the viral parameter. This avoids a 
continuously drifting star formation rate such as in the KM 
simulation. Nonetheless, the level of e mo d following from the 
PHN model is systematically too low in the regime that is 
relevant for star formation. As we will show in Section [3.3| 
this results in an extremely clumpy disk, which only slowly 
reaches a stable configuration. Although the fit formula for 
£mod put forward by prof] depends only on the virial pa¬ 
rameter <a v ir oc Ksgs, the lower envelope of the e mo d-A"sGS- 
distribution shown in Fig.[2]suggests a quadratic dependence 
on i^sGS- This non-linear behaviour is a consequence of the 
implicit iFsGS-dependencies of the cold-gas density p c , P a and 
the clump scale £ c . 

3.3 Disk structure and probability density 
functions 

By comparing the gaseous disk structure for the different 
runs in Fig. [ 3 ] it can be seen that deviations of the star 
formation rate caused by different average levels of e mo d go 
along with changes in the disk structure. This is obvious 
for the cases KM, KM-FK, and PHN, which produce much 
clumpier disks. To a lesser degree, also the disk in the PN 
run tends to be clumpier as compared to PN-FK. Virtually 
no differences can be seen for the runs from the HC-family, 
which also compare well to PN-FK. In extreme cases (KM 
and PHN), the disks are dominated by the presence of very 
massive, and dense clumps (E > 10 3 M 0 pc -2 ). The void 
areas (S < 1 M 0 pc -2 ) between these clumps are only in¬ 
tercepted by tidal tails and bridges that form during mergers 
and fly-by interactions of clusters. As a result, a similar star 
formation rate such as in the late stages of the PHN and 
PN-FK runs does not necessarily imply similar disk struc¬ 
tures. 

In the KM and PHN runs, we find extremely massive 

5 This formula is not based on any physical ideas. It is merely a 
fit to the data points obtained from their simulation suite. 


clumps with very high column densities (E > 10 3 M 0 pc -2 ). 
This is in stark contrast to other runs (PN-FK, PN-FK-mff, 
HC family, and KM-FK-mff), where gas of moderate density 
(1 M 0 pc -2 < E < 100 M 0 pc -2 ) is concentrated in spiral¬ 
like structures with substantially smaller clumps, similar to 
the initial disk at 1 Gyr. Denser, star-forming regions are 
eventually disrupted by the violent expansion induced by 
SN feedback, but the low-density gas is less space filling in 
these runs. The disks in the PN and the KM-FK runs are 
intermediate cases. 

A more accurate distinction is possible by means of 
probability density functions of the gas column density 
weighted by by area or by the mass fraction Ex/E, where 
Ex either stands for the column density of cold gas, E c , or 
the column density of shielded molecular gas, Eh 2 - The re¬ 
sults for the different simulation runs are shown in Fig. 
The normalisations of pdf [log 10 (E)] and pdf Ex / E [log 10 (E)] 
are, respectively, given by 

/ oo 

pdf[log 10 (S)]d[log 10 (S)] = 1 (14) 

-oo 

and 

pdf [log 10 (E)]d[log 10 (S)] = (15) 

One can distinguish at least four different regimes in the 
plot of pdf[log 10 (E)] in the upper left panel of Fig. [4] 

(i) E < 2 x 10 -1 M 0 pc -2 : Gas found in this range is 
extremely hot and dilute and was most likely recently pro¬ 
duced by SNe. 

(ii) 2 x 10 -1 M 0 pc -2 < E < 3 x 10° M 0 pc -2 : In this 
range and also at higher densities, there is an approxi¬ 
mate balance between heating and radiative cooling. The 
two prominent maxima of the pdfs correspond to the col- 
lisional ionisation of hydrogen at ~ 2 X 10° M 0 pc -2 and 
helium at ~ 5 x 10 -1 M 0 pc -2 , respectively. For this reason, 
there is hardly any gas in the cold phase below a density of 
~ 2 x 10° M 0 pc -2 (compare top and middle plots in Fig. [ 4 ]). 

(iii) 3 x 10° M 0 pc -2 < E < 10 2 M 0 pc -2 : In this 
regime, the cold and warm phases coexist. The distribu¬ 
tion of the cold gas shown in the middle plot in Fig. [4] 
has an approximately log-normal shape, with a maximum 
at E « 2 x 10 1 M 0 pc -2 . The existence of a cold phase can 
also be seen as a hump in pdf[log 10 (E)]. Above densities of 
about « 1 x 10 1 M 0 pc -2 , shielded molecular gas may be 
present inside of the cold phase (see bottom plot). 

(iv) E > 10 2 M 0 pc -2 : Due to significant gravitational 
self-interactions of gas and stars, both pdf [log 10 (E)] and 
pdf Sc / E [log 10 (E)] show an excess in their high density tails. 
The gas in this regime resides in the prominent, more or less 
gravitationally bound knots and clumps that are visible in 
Fig. Shielded molecular gas makes up the major fraction 
of the cold phase, but there is also gas in the warm phase 
because of the evaporation of cold gas due to SN feedback. 

Differences induced by the star formation models are 
most evident in the density regime (iv), in which most of 
the star formation happens. There is an enormous excess 
of high-density gas in case of PHN and KM. The maxima 
of shielded molecular gas are found at roughly ten times 
higher E (« 6 X 10 1 M 0 pc -2 ) compared to the maxima in 
the other runs. The lack of intermediate density gas as well 
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Figure 3. Projections of the total gas density p perpendicular 
to the disk plane in an area of 15 X 15 kpc around the galactic 
centre. The individual projections depict the state of the disk in 
the indicated run at ~1.4 Gyr. 


as the increased abundance of low-density gas in PHN and 
KM is readily visible in the regimes iii) and iv). The pdfs 
from the KM-FK and the PN run exhibit significantly less 
pronounced high-density tails, but they are still stronger 
than the tails obtained with the HC-family, PN-FK, PN- 
FK-mff, and KM-FK-mff models. 

Observationally almost log-normal shaped pdfs are 


found for gas in nearby galaxies (e.g. Berkhuijsen & Fletcher 
2015 Hughes et al.||2013| . The pdf of CO-line emitting gas 


in M51 by Hughes et al. ( |2013| is based on data of compa¬ 
rable spatial resolution (40 pc). Most of what we consider 
as cold phase gas (E c ) would be seen in such observations, 
but a significant fraction of it is actually atomic, particu¬ 
larly for lower E. Compared to observed pdfs (e.g. Figure 2a 


in Hughes et al. 2013), the maximum of pdf s ,/£ l°glo( S )] 
computed from our simulations is thus shifted toward lower 
densities. Apart from that, M51 is an interacting galaxy in 
which highly compressed gas is expected to be more abun¬ 
dant than in an isolated, quiescent galaxy. Nevertheless, the 
general shape and especially the high-density tail as inferred 
from the HC-family, PN-FK, PN-FK-mff, and KM-FK-mff 
model runs are consistent with the observations of |Hughes| 


et al. (2013), while PHN and KM are clearly ruled out. 


A comparison with the simulations of |Hopkins et al.| 


(2012) shows pdfs that are populated up to much higher 
densities. This appears to be a consequence of the signifi¬ 
cantly higher numerical resolution, but also the star forma¬ 
tion recipe used in their simulations. Specifically, the low ef- 



£ [M 0 /pc 2 ] £ [M 0 /pc 2 ] 


w o 

1 



KT 1 



£[M e /pc 2 ] £[M 0 /pc 2 ] 


■ ■ ■ PN ■ - HC - - KM - - PHN 

- PN-FK - HC-FK - KM-FK 

- - PN-FK-mff - - HC-FK-mff - - KM-FK-mff 


Figure 4. Top left: Area weighted probability density func¬ 
tion pdf(log 10 E) of the logarithmic total gas column density in 
the disks versus total gas column density S for all runs at a 
simulated age of ~ 1.4 Ga. Arrangement of line colours and 
styles like in Fig. [T] Top right: Detail of the plot on the left 
for E > 10 2 M©pc -3 . Middle left: Probability density func¬ 
tion pdf Sc/ / E (log 10 S) weighted by the cold phase gas fraction 
E c /E versus E. Middle right: Detail of the plot on the left for 
S > 10 2 Mqpc -3 . Bottom left: Probability density function 
pdf Sc/ / E (log 10 E) weighted by the shielded molecular gas frac¬ 
tion Eh 2 /E versus E. Bottom right: Detail of the plot on the 
left for E > 10 2 Mqpc -3 . 


ficiency and the rather high star formation threshold favour 
the formation of denser clouds and a higher global star for¬ 
mation rate (cf. their ’Sbc’-runs). 


3.4 Star formation relations 

With the exception of the KM and PHN runs, a robust 
relation between the star-format ion column density Esf 
and column density of shielded, molecular gas Eh 2 is es¬ 
tablished (see Fig. [5]). However, the inferred time scales 
tsf,h 2 of star formation in molecular gas vary between 
« 80 Myr and « 0.5 Gyr. As pointed out in |BS14[ time 
scales around ~ 100 Myr are reasonable because the simpli¬ 
fied model for the molecular fraction /h 2 in our simulations 
is biased toward high-density, star-forming regions. Observa¬ 
tional analogues of what is treated as shielded molecular gas 
in our simulations favour short star formation time scales 
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molecular surface density S H [M 0 pc 2 ] 



Figure 5. Comparison of area-weighted 2D histograms of the 
star formation column density E s vs. the column density Eh 2 
of shielded, molecular gas at an age of ~ 1.4 Gyr. The area is 
logarithmically colour-coded. 


100 Myr (e.g. Gao &; Solomon 2004 Evans et al. 2009 
Murray 2011). Keeping that in mind, the star formation rate 


is anomalously low for a wide range of Eh 2 column densi¬ 
ties for the KM- and the PHN-model, with a very sharp rise 
at Eh 2 ~ 10 3 M©pc -2 . Moreover the spread is significantly 
larger in these cases. Such Esf-Eh 2 -relations are definitely 
not consistent with observational counterparts. For the KM- 
FK and PN models, the relations are slightly tilted upwards 
at the high-density end|^] Basically, the results show that 
star formation is shifted to higher molecular column densi¬ 
ties for the models in the KM family, particularly for the 
original KM model, and also for the PHN model. Gener¬ 
ally, the multi free-fall models produce tighter correlations 
between Esf and Eh 2 . The tightest correlation is found for 
HC-FK-mff, with a particularly short star formation time 
scale. 

The 2D histograms of Esf vs. E (total gas column 
density) in Fig. |6| sh ow qualitatively the same behaviour 
as described in IBS14I Star formation is cut off below E « 
10° M©pc -2 because of the lack of shielded, molecular gas. 


6 This is actually a consequence of a lower bound for the cold-gas 
fraction at densities above 500 cm -3 . For the fiducial model, this 
was intended to prevent star formation from being completely 
quenched by feedback in very dense regions, which hardly af¬ 
fected the Esf~Eh 2 -relation (see for example, the PN-FK panel 
in Fig. [5j. For extremely clumpy disks, the star formation rate 
is artificially pushed upwards. We did not adjust the threshold 
density for the cold-gas floor when carrying out our comparison 
study. 


Figure 6. 2D histograms as in Fig. [5] for the star formation col¬ 
umn density E s vs. the total gas column density E. 


Above this threshold, Esf steeply rises and then gradually 
flattens to an approximate KS-relation Esf oc Eg s with 
c^ks ~ 2 . An exponent ckks ~ 2 is indicative of star bursts. 
This is plausible, since the simulated galaxies are still gas- 
dominated after 1.4 Gyr. Again, EsF-E-relations obtained 
with the models KM, PHN, and (to a lesser degree) KM- 
FK differ for the same reasons like in the context of the 
Esf-Eh 2 -relations. 


4 DISCUSSION AND CONCLUSIONS 


In this paper we presented results of a suite of ten isolated 
disk galaxy simulations carried out using the cosmological 
hydrodynamics code Nyx (Alm gren et al.||2013 ), featuring a 
turbulence SGS model ( Schmidt fe Federrath|20lT Schmidt 
et al.|20l4 ) and the MIST model (BS12 BS14) for the tur¬ 
bulent multi-phase ISM. In each of the runs, a different an¬ 
alytic star-formation model from the literature was applied 
flFK12| |HC08-13| |KM05| |PN11| |PHN12[ see Table [l] for the 
analytical expressions employed here). These models for the 
sub-resolution structure of the ISM were incorporated into 
the MIST framework to estimate the local star formation ef¬ 
ficiency e mo d in the cold molecular gas. We used an evolved 
(for 1 Gyr) disk from the ’ref’ run of BS14[ which resem¬ 
bles a gas-rich, quiescent star-forming spiral galaxy, as ini¬ 
tial conditions to start the runs with different star forma¬ 
tion models. This approach minimised effects caused by the 
rapid growth of stellar mass and metallicity during the ini¬ 
tial transient. From these initial conditions every simulation 
was evolved over additional 0.4 Gyr, which allows the star 
forming disk to settle into a new equilibrium in most cases. 
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10 H. Braun and W. Schmidt 


Several important properties resulting from the local self¬ 
regulation in MIST, as reported by |BS14[ are recovered in 
our simulation runs independent of the choice of the star for¬ 
mation efficiency e mo d if the KM, KM-FK, and PHN models 
are excluded: 

• A global star formation rate around Msf ~ 2.5 M©yr -1 
is maintained. This is a typical value for a quiescent isolated 
spiral galaxy. 

• Model-dependent variations of e mo d do not significantly 
affect the amount of shielded molecular gas. As a result, 
the almost linear relationship between local star formation 
rate column density E s and column density Eh 2 of shielded 
molecular gas turns out to be very robust in the frame¬ 
work of MIST. However, the inferred depletion time scales 
vary between 80 Myr and 0.5 Gyr, depending on the star 
formation model. With respect to observational findings, a 


depletion time scale < 100 Myr is more plausible (Gao & 
Solomon 2004| Murray 2011). Moreover, a tighter correla¬ 
tion is obtained for the multi free-fall models. 

• Our star-forming regions have a typical life-cycle of 
about 10 to 30 Myr. A fraction of ~ 10 per cent of their 
gaseous mass is converted into stars during that period, 


which matches observational estimates (e.g. Blitz et al. 2007 


McKee &; Qstriker|2007 Miura et al.|2012 ) as well as simu¬ 

lation results by e.g. Hopkins et al. ( |2012| ) and Tasker et al. 
(2015). We also find long-lived giant star-forming clumps 
similar to those described by Bournaud et al. (2014) in our 
simulations. 

• Cold-gas density pdfs are comparable to observational 
counterparts (see, for example, Hughes et al. 2013]), partic¬ 
ularly at the high densities that are required for star forma¬ 
tion. 

• The high-density tails of our E s -E-distributions suggest 
a KS-relationship E s oc E as with ~ 2, which is expected 
for environments dominated by gas dynamics. 

• Star formation mainly occurs in the SGS turbulence 
energy range 10 km 2 s' 2 < Ksgs < 10 3 km 2 s 2 , with a 
peak around 10 2 km 2 s -2 . This value corresponds to a 3D- 
velocity dispersion in the range between cr u ~ 3 kms -1 and 
30 kms -1 , which is in good agreement with observations of 
star-forming regions (e.g. | Leroy etT~a l. 2008 Sh etty et al. 
2012 Stilp et ak]|2013| (since these regions are more or less 
unresolved in our simulations, the SGS turbulence energy 
T^sgs is the relevant quantity). A minimum level of turbu¬ 
lence is necessary to boost the molecular hydrogen formation 
rates. On the other side, there is an upper bound on Ksgs 
because high values of Ksgs, which are produced by intense 
stellar feedback, tend to disrupt star forming regions. This 
effect also limits the duty cycle of star-forming regions to 
around 10 to 30 Myr in our simulation runs, in agreement 
with observational estimates (e.g. Blitz et al.||2007 McKee 


&; Qstriker|2007| Miura et al.|2012). 


For most models, e mo d is effectively a function of Ksgs- The 
scatter of e mo d in Fig. [2] indicates a subdominant depen¬ 
dency on other quantities such as the density, for example, 
in the case of PHN. For e mo d ~ 0.1, corresponding to molec¬ 
ular gas depletion time scales of < 100 Myr in the simula¬ 
tions, the gaseous star forming disks are moderately clumpy 
and exhibit pronounced transient spiral-like, flocculent fea¬ 
tures connecting knots and clumps. If e mo d is substantially 
lower in the interval of turbulence energy values that admit 


significant star formation, the disk structure changes. This 
particularly happens with the KM, KM-FK, and PHN mod¬ 
els, which produce either too many or not enough stars to 
support the initial disk configuration through feedback. In 
these runs, the disk undergoes a transition to a different con¬ 
figuration with markedly more massive and compact clumps 
accompanied by tidal tails and bridges. Such extremely mas¬ 
sive clumps appear to be common in star-forming, gas-rich 

Zanella et al.||2015 Guo 


galaxies at high redshifts 


et al.|2012 ). Similar structures are found in some simulations 
of gas-rich galaxies (e.g. by Hopkins et al. 2011 Bournaud 


et al. 2014). However, the models considered here are in¬ 


tended to describe conditions in fairly regular star-forming 
regions, which are found in almost bulge-less and gas-rich 
spiral galaxies without any external disturber. Some of the 
models clearly fail to produce a disk structure that even 
remotely resembles such galaxies. However, the strong neg¬ 
ative coupling between star formation efficiency and turbu¬ 
lence produced by stellar feedback, for example, in the case 
of the KM model, might help to gain a better understand¬ 
ing of galaxies dominated by massive clumps. Which factors 
actually select the mode of star formation is left as an open 
question. 

If additional feedback processes like winds from young, 
massive stars were incorporated in our simulations, the 
global star formation rate would be shifted, but the self¬ 
regulation mechanism would basically remain intact. Apart 
from that, the star formation rate would be altered if the cir- 
cumgalactic medium were properly taken into account (such 
as in cosmological zoom-in simulations). This has to be in¬ 
vestigated in future studies. 

Both KM and KM-FK do not feature a free-fall time 
factor (i.e. /s = 1, meaning that the effective free-fall time 
scale is given by the mean density of the gas in cold phase), 
while all other theoretical models use some weighting fac¬ 
tor to modify the effective free-fall time scale. In our sim¬ 
ulations, this turns out to be an essential feature to repro¬ 
duce observed properties of nearby disk galaxies and their 
star-forming regions. The importance of a non-trivial free- 
fall factor, such as in the PN and the HC models, was also 
highlighted by |FK12[ Moreover, we can confirm that the 
calibrations established by this study generally result in im¬ 


provements. Although 


FK12 


found higher residual x -values 


for the best-fit HC models than for the competing models, 
our simulations suggest that the models from the HC family 
yield the most plausible results in terms of the correlation 
between star formation rate and molecular hydrogen column 
density, the cold-gas pdf, and the disk structure. However, 
the differences compared to the PN-FK and PN-FK-mff runs 
are minor. Also the KM model produces consistent results if 
it is applied in the multi free-fall formulation with the new 
calibration of ilFK12 

In contrast to all other models in our suite, the star 
formation efficiency in the simple PHN model has no func¬ 
tional dependence on the Mach number and, consequently, 
is ignorant of the fundamentally different nature of subsonic 
and supersonic turbulent fluctuations. This does not appear 
plausible, because the width of the density pdf is basically 
determined by the turbulent Mach-number in the cold phase 
M (see equation [6|. However, PHN12 incorporate magnetic 
fields, which might also influence the star formation effi¬ 
ciency, although probably to a lesser degree. Possibly, the 
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model could be improved by establishing a fit function from 
turbulence simulations covering a larger parameter space, 
e.g. by including non-solenoidal forcing. 

A further distinction between the models that are so far 
consistent with the observational constraints - namely, PN- 
FK, PN-FK-mff, HC, HC-FK, HC-FK-mff, and KM-FK-mff 
- might be possible on the basis of observational measure¬ 
ments of the turbulent velocity dispersion on molecular- 
cloud scales and associated star formation efficiencies. On 
the theoretical side, constraints on the relation between star 
formation efficiency and parameters related to turbulence 
can be obtained from high-resolution simulations of star¬ 
forming clouds such as in |FK12] However, any such simula¬ 
tions make idealised assumptions about initial and bound¬ 
ary conditions for molecular-cloud turbulence. To check the 
consistency of models for the dynamical star formation ef¬ 
ficiency, it will be necessary to re-simulate a star-forming 
region in global disk galaxy simulation such that star for¬ 
mation processes are fully resolved and then to compare the 
results with the assumptions that went into the model. 
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